# Mark Nieman and Maxwell Allamong
# School of Thought - Data Analysis
# November 7, 2022
sink("NiemanAllamong_SoT.txt")
rm(list = ls())

#Load Packages
library(tidyverse)
library(magrittr)
library(stargazer)
library(sandwich)


# Load data ----
SoT_data<-read_csv("SoT_data.csv") 

###################
## DATA ANALYSIS ##
###################


# Table 1: Leader Educational Background
Tab1 <- SoT_data 
Tab1$AngloAmer.edu[Tab1$uni.type == "Military" & (Tab1$uni.countryabb == "USA" | Tab1$uni.countryabb == "UKG" | Tab1$uni.countryabb == "AUL" | Tab1$uni.countryabb == "NEW" | Tab1$uni.countryabb == "CAN")]<-1
Tab1$euro.edu[Tab1$uni.type == "Military" & (Tab1$uni.countryabb == "AUS" | Tab1$uni.countryabb == "BEL" | Tab1$uni.countryabb == "DEN" |
                                               Tab1$uni.countryabb == "FIN" | Tab1$uni.countryabb == "FRN" | 
                                               Tab1$uni.countryabb == "GMY" | Tab1$uni.countryabb == "GRC" | 
                                               Tab1$uni.countryabb == "IRE" | Tab1$uni.countryabb == "ITA" | 
                                               Tab1$uni.countryabb == "NOR" | Tab1$uni.countryabb == "NTH" | 
                                               Tab1$uni.countryabb == "POR" | Tab1$uni.countryabb == "SPN" | 
                                               Tab1$uni.countryabb == "SWD" | Tab1$uni.countryabb == "SWZ")] <-1
Tab1$noneuro.edu[Tab1$uni.type == "Military" & Tab1$AngloAmer.edu==0 & Tab1$euro.edu==0]<-1
Tab1$education[Tab1$none.edu==1]<-"None"
Tab1$education[Tab1$AngloAmer.edu==1]<-"Anglo-American"
Tab1$education[Tab1$euro.edu==1]<-"Other Western"
Tab1$education[Tab1$noneuro.edu==1]<-"Hierarchical"
Tab1$abroad.edu<-NA
Tab1$abroad.edu[Tab1$uni.countryabb!=Tab1$idacr & (Tab1$mil.edu | Tab1$education!="None")]<-"Abroad"
Tab1$abroad.edu[Tab1$uni.countryabb==Tab1$idacr & (Tab1$mil.edu | Tab1$education!="None")]<-"Domestic"
table(Tab1$education,Tab1$abroad.edu,Tab1$mil.edu)


#Table 2: Regional Distribution of Leader Education Background
Tab1$aa[Tab1$uni.countryabb=="USA"]<-"USA"
Tab1$aa[Tab1$uni.countryabb=="UKG"]<-"UK"
Tab1$aa[Tab1$uni.countryabb=="NEW"]<-"New Zealand"
Tab1$aa[Tab1$uni.countryabb=="AUL"]<-"Australia"
Tab1$aa[Tab1$uni.countryabb=="CAN"]<-"Canada"
Tab1$aa[Tab1$euro.edu==1 ]<-"Other Western"
Tab1$aa[Tab1$noneuro.edu==1]<-"Hierarchical"
Tab1$aa[Tab1$none.edu==1]<-"None"
table(Tab1$aa)
table( Tab1$aa, Tab1$region)
mil <- Tab1 %>% filter(uni.type=="Military")
table( mil$aa,mil$region)

#Table 3: Leader Subject Studied at University
table(SoT_data$AA.econ)
table(SoT_data$AA.law)
table(SoT_data$AA.socsci)
table(SoT_data$AA.hum)
table(SoT_data$AA.other)
table(SoT_data$noAA.econ)
table(SoT_data$noAA.law)
table(SoT_data$noAA.socsci)
table(SoT_data$noAA.hum)
table(SoT_data$noAA.other)


# Main Analysis: Liberalization Models for Tables 4 and 5

# Trade Liberalization
# Tertiary Education 
SoT_data1<- SoT_data %>% filter(start.trade.liberal == 0)
trade1 <- glm(change.trade.liberal ~ AngloAmer.edu + euro.edu + noneuro.edu  +
                           log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure, data = SoT_data1, family = "binomial", na.action=na.exclude)
trade1.se <- vcovCL(trade1, type = "HC0", SoT_data1$ccode)
  # Wald tests
  (coef(trade1)[2]- coef(trade1)[3])/sqrt(trade1.se[2,2]+trade1.se[3,3]-2*trade1.se[2,3])
  (coef(trade1)[2]- coef(trade1)[4])/sqrt(trade1.se[2,2]+trade1.se[4,4]-2*trade1.se[2,4])
trade1.se <- sqrt(diag(trade1.se))

trade1fe <- glm(change.trade.liberal ~ AngloAmer.edu + euro.edu + noneuro.edu +
                          factor(decade) + factor(region), data = SoT_data1, family = "binomial")
trade1fe.se <- vcovCL(trade1fe, type = "HC0", SoT_data1$ccode)
  # Wald tests
  (coef(trade1fe)[2]- coef(trade1fe)[3])/sqrt(trade1fe.se[2,2]+trade1fe.se[3,3]-2*trade1fe.se[2,3])
  (coef(trade1fe)[2]- coef(trade1fe)[4])/sqrt(trade1fe.se[2,2]+trade1fe.se[4,4]-2*trade1fe.se[2,4])
trade1fe.se <- sqrt(diag(trade1fe.se))

# Specialization
trade2 <- glm(change.trade.liberal ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu +
                             log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                             EthFrac + log(pop) + tenure, data = SoT_data1, family = "binomial", na.action=na.exclude)
trade2.se <- vcovCL(trade2, type = "HC0", SoT_data1$ccode)
trade2.se <- sqrt(diag(trade2.se))

trade2fe <- glm(change.trade.liberal ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu +
                               factor(decade) + factor(region), data = SoT_data1, family = "binomial")
trade2fe.se <- vcovCL(trade2fe, type = "HC0", SoT_data1$ccode)
trade2fe.se <- sqrt(diag(trade2fe.se))
  

# Rule of Law 
# Tertiary Education
rl1 <- lm(change.LJI ~ AngloAmer.edu + euro.edu + noneuro.edu +
                           log(rgdppc) + oil +  UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.LJI, data = SoT_data, na.action=na.exclude)
rl1.se <- vcovCL(rl1, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(rl1)[2]- coef(rl1)[4])/sqrt(rl1.se[2,2]+rl1.se[4,4]-2*rl1.se[2,4])
rl1.se <- sqrt(diag(rl1.se))
  
rl1fe <- lm(change.LJI ~ AngloAmer.edu + euro.edu + noneuro.edu + start.LJI + 
                     factor(decade) + factor(region), data = SoT_data)
rl1fe.se <- vcovCL(rl1fe, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(rl1fe)[2]- coef(rl1fe)[3])/sqrt(rl1fe.se[2,2]+rl1fe.se[3,3]-2*rl1fe.se[2,3])
  (coef(rl1fe)[2]- coef(rl1fe)[4])/sqrt(rl1fe.se[2,2]+rl1fe.se[4,4]-2*rl1fe.se[2,4])
rl1fe.se <- sqrt(diag(rl1fe.se))
  
# Specialization 
rl2 <- lm(change.LJI ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu +
                           log(rgdppc) + oil +  UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.LJI, data = SoT_data, na.action = na.exclude)
rl2.se <- vcovCL(rl2, type = "HC1", SoT_data$ccode)
rl2.se <- sqrt(diag(rl2.se))

rl2fe <- lm(change.LJI ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu + start.LJI + 
                             factor(decade) + factor(region), data = SoT_data)
rl2fe.se <- vcovCL(rl2fe, type = "HC1", SoT_data$ccode)
rl2fe.se <- sqrt(diag(rl2fe.se))
  

# Financial liberalization 
# Tertiary Education
ka1 <- lm(change.kaopen ~ AngloAmer.edu + euro.edu + noneuro.edu + 
                           log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.kaopen, data = SoT_data, na.action=na.exclude)
ka1.se <- vcovCL(ka1, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(ka1)[2]- coef(ka1)[3])/sqrt(ka1.se[2,2]+ka1.se[3,3]-2*ka1.se[2,3])
  (coef(ka1)[2]- coef(ka1)[4])/sqrt(ka1.se[2,2]+ka1.se[4,4]-2*ka1.se[2,4])
  (coef(ka1)[3]- coef(ka1)[4])/sqrt(ka1.se[3,3]+ka1.se[4,4]-2*ka1.se[3,4])
ka1.se <- sqrt(diag(ka1.se))

ka1fe <- lm(change.kaopen ~ AngloAmer.edu + euro.edu + noneuro.edu + start.kaopen + 
                     factor(decade) + factor(region), data = SoT_data)
ka1fe.se <- vcovCL(ka1fe, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(ka1fe)[2]- coef(ka1fe)[3])/sqrt(ka1fe.se[2,2]+ka1fe.se[3,3]-2*ka1fe.se[2,3])
  (coef(ka1fe)[2]- coef(ka1fe)[4])/sqrt(ka1fe.se[2,2]+ka1fe.se[4,4]-2*ka1fe.se[2,4])
  (coef(ka1fe)[3]- coef(ka1fe)[4])/sqrt(ka1fe.se[3,3]+ka1fe.se[4,4]-2*ka1fe.se[3,4])
ka1fe.se <- sqrt(diag(ka1fe.se))
  
# Specialization
ka2 <- lm(change.kaopen ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu +
                           log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.kaopen, data = SoT_data, na.action=na.exclude)
ka2.se <- vcovCL(ka2, type = "HC1", SoT_data$ccode)
ka2.se <- sqrt(diag(ka2.se))
  
ka2fe <- lm(change.kaopen ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu + start.kaopen + 
                             factor(decade) + factor(region), data = SoT_data)
ka2fe.se <- vcovCL(ka2fe, type = "HC1", SoT_data$ccode)
ka2fe.se <- sqrt(diag(ka2fe.se))
  

# Human Rights
# Tertiary Education
hr1 <- lm(change.hrscore ~ AngloAmer.edu + euro.edu + noneuro.edu +
                           log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.hrscore, data = SoT_data, na.action=na.exclude)
hr1.se <- vcovCL(hr1, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(hr1)[2]- coef(hr1)[4])/sqrt(hr1.se[2,2]+hr1.se[4,4]-2*hr1.se[2,4])
hr1.se <- sqrt(diag(hr1.se))

hr1fe <- lm(change.hrscore ~ AngloAmer.edu + euro.edu + noneuro.edu + start.hrscore+ 
                            factor(decade) + factor(region), data = SoT_data)
hr1fe.se <- vcovCL(hr1fe, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(hr1fe)[2]- coef(hr1fe)[4])/sqrt(hr1fe.se[2,2]+hr1fe.se[4,4]-2*hr1fe.se[2,4])  
hr1fe.se <- sqrt(diag(hr1fe.se))
  
# Specialization
hr2 <- lm(change.hrscore ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu + 
                           log(rgdppc) + oil + xconstV + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.hrscore, data = SoT_data, na.action=na.exclude)
hr2.se <- vcovCL(hr2, type = "HC1", SoT_data$ccode)
hr2.se <- sqrt(diag(hr2.se))
  
hr2fe <- lm(change.hrscore ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu + start.hrscore  + 
                             factor(decade) + factor(region), data = SoT_data)
hr2fe.se <- vcovCL(hr2fe, type = "HC1", SoT_data$ccode)
hr2fe.se <- sqrt(diag(hr2fe.se))
  

# Liberal Democracy
# Tertiary Education
ld1 <- lm(change.libdem ~ AngloAmer.edu + euro.edu + noneuro.edu + 
                           log(rgdppc) + oil +  UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.libdem, data = SoT_data, na.action=na.exclude)
ld1.se <- vcovCL(ld1, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(ld1)[2]- coef(ld1)[3])/sqrt(ld1.se[2,2]+ld1.se[3,3]-2*ld1.se[2,3])
  (coef(ld1)[3]- coef(ld1)[4])/sqrt(ld1.se[3,3]+ld1.se[4,4]-2*ld1.se[3,4])    
ld1.se <- sqrt(diag(ld1.se))
  
ld1fe <- lm(change.libdem ~ AngloAmer.edu + euro.edu + noneuro.edu + start.libdem + 
                             factor(decade) + factor(region), data = SoT_data)
ld1fe.se <- vcovCL(ld1fe, type = "HC1", SoT_data$ccode)
  # Wald tests
  (coef(ld1fe)[2]- coef(ld1fe)[4])/sqrt(ld1fe.se[2,2]+ld1fe.se[4,4]-2*ld1fe.se[2,4])   
  (coef(ld1fe)[3]- coef(ld1fe)[4])/sqrt(ld1fe.se[3,3]+ld1fe.se[4,4]-2*ld1fe.se[3,4])   
ld1fe.se <- sqrt(diag(ld1fe.se))
  
# Specialization
ld2 <- lm(change.libdem ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu +
                           log(rgdppc) + oil  + UK_col + us_ally + 
                           EthFrac + log(pop) + tenure + start.libdem, data = SoT_data, na.action=na.exclude)
ld2.se <- vcovCL(ld2, type = "HC1", SoT_data$ccode)
ld2.se <- sqrt(diag(ld2.se))
  
ld2fe <- lm(change.libdem ~ AA.econ + AA.law + AA.socsci + AA.hum + AA.other + noAA.econ + noAA.law+ noAA.socsci + noAA.hum + noAA.other + AAmil.edu + noAAmil.edu + start.libdem + 
                             factor(decade) + factor(region), data = SoT_data)
ld2fe.se <- vcovCL(ld2fe, type = "HC1", SoT_data$ccode)
ld2fe.se <- sqrt(diag(ld2fe.se))
  

# Table 4: Tertiary Educational Model and Political Outcomes
stargazer(trade1, trade1fe, rl1, rl1fe, ka1, ka1fe, hr1, hr1fe, ld1, ld1fe, digits = 3,
            se = list(trade1.se, trade1fe.se, rl1.se, rl1fe.se, ka1.se, ka1fe.se, hr1.se, hr1fe.se, ld1.se, ld1fe.se),
            intercept.bottom = T, dep.var.labels = c("Trade Liberalization","Rule of Law","Financial Openness","Human Rights","Liberal Democracy"),
            covariate.labels = c("Anglo-American Education", "Other Western Education", "Hierarchical Education", 
            "Economic Development","Oil Producer","Executive Constraints","Former British Colony", "US Ally",
            "Ethnic Fractionalization","Population","Time in Office","Entered 1960s","Entered 1970s","Entered 1980s","Entered 1990s","Entered 2000s", "Entered 2010s",
            "Asia", "Eastern Europe","MENA","Oceania", "Subsaharan Africa", "LJI Entering Office", "Entered 1950s", "KA Entering Office", "HR Entering Office", "LD Entering Office"),
            column.labels = c("Edu. Model", "FEs", "Edu. Model", "FEs", "Edu. Model", "FEs", "Edu. Model", "FEs", "Edu. Model", "Fixed Effects"),
            keep.stat = c("n","ll","rsq"),align=T,no.space = T,
            star.char = c("+","*", "**"),
            star.cutoffs = c(0.2, 0.1, 0.05),         
            notes = c("$^{**}$p$<$0.05$, ^{*}$p$<$0.1, two-tailed; $^{+}$p$<$0.1, one-tailed;. Standard errors clustered by country."), 
            notes.append = F)
  
# Table 5: Tertiary Educational Model, Subject Studied, and Political Outcomes
stargazer(trade2, trade2fe, rl2, rl2fe, ka2, ka2fe, hr2, hr2fe, ld2, ld2fe, digits = 3,
            se = list(trade2.se, trade2fe.se, rl2.se, rl2fe.se, ka2.se, ka2fe.se, hr2.se, hr2fe.se, ld2.se, ld2fe.se),
            intercept.bottom = T, dep.var.labels = c("Trade Liberalization","Rule of Law","Financial Openness","Human Rights","Liberal Democracy"),
            covariate.labels = c("AA Economics", "AA Law", "AA Social Science", "AA Humanities", "AA Other Majors",   
                                 "Non-AA Economics", "Non-AA Law", "Non-AA Social Science", "Non-AA Humanities", "Non-AA Other Majors",
                                 "AA Military School","Non-AA Military School",
                                 "Economic Development","Oil Producer","Executive Constraints","Former British Colony", "US Ally",
                                 "Ethnic Fractionalization","Population","Time in Office",
                                 "Entered 1960s","Entered 1970s","Entered 1980s","Entered 1990s","Entered 2000s", "Entered 2010s",
                                 "Asia", "Eastern Europe","MENA","Oceania", "Subsaharan Africa", "LJI Entering Office", "Entered 1950s", "KA Entering Office", "HR Entering Office", "LD Entering Office"),
            column.labels = c("Subject", "FEs", "Subject", "FEs", "Subject", "FEs", "Subject", "FEs", "Subject", "Fixed Effects"),
            keep.stat = c("n","ll","rsq"),align=T,no.space = T,
            star.char = c("+","*", "**"),
            star.cutoffs = c(0.2, 0.1, 0.05),         
            notes = c("$^{**}$p$<$0.05$, ^{*}$p$<$0.1, two-tailed; $^{+}$p$<$0.1, one-tailed;. Standard errors clustered by country."), 
            notes.append = F)
  

# Figure 1: Two Dimensions of Tertiary Education Models
library(ggplot2)
mylabs<-c("Soviet","Post-Soviet","East Asian","Latin","Humboldtian","Bonapartist","Anglo-American")
hierarchy<-c(4,4,4,2,1.5,2,0)
state<-c(1,3.5,1.5,2,5.5,4.5,7)

ggplot()+geom_text(aes(x=(hierarchy)+1, y=state, label=mylabs), size=5)+
  geom_vline(aes(xintercept=3),linetype=2, color="gray")+
  geom_hline(aes(yintercept=4),linetype=2, color="gray")+
  scale_x_continuous("Hierarchy in the Classroom",limits = c(0,6), 
                     breaks=c(1.3,4.8),labels=c("Egalitarian","Authoritarian"),position = "top")+
  scale_y_continuous("Autonomy from the State",breaks=c(2.25,6.25),limits = c(0.5,7.5),labels=c("Low","High"))+
  theme_bw()+
  theme(axis.title = element_text(color="black", size=16, face="bold"),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size=12,angle=90, face="bold"),
        axis.text.x = element_text(size=12, face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

ggsave(file = "Uni_types.pdf", height=6, width=8)

